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Abstract 

As nuclear wave functions have to obey the Pauli principle, potentials issued from reaction theory or Hartree- 
Fock formalism using finite-range interactions contain a non-local part. Written in coordinate space represen- 
tation, the Schrodinger equation becomes integro-differential, which is difficult to solve, contrary to the case of 
local potentials, where it is an ordinary differential equation. A simple and powerful method has been proposed 
several years ago, with the trivially equivalent potential method, where non-local potential is replaced by an 
equivalent local potential, which is state-dependent and has to be determined iteratively. Its main disadvantage, 
however, is the appearance of divergences in potentials if the wave functions have nodes, which is generally the 
case. We will show that divergences can be removed by a slight modification of the trivially equivalent potential 
method, leading to a very simple, stable and precise numerical technique to deal with non-local potentials. 
Examples will be provided with the calculation of the Hartree-Fock potential and associated wave functions of 
^^O using the finite-range N^LO realistic interaction. 

PACS 

02.60.Nm Integral and integrodifferential equations 
03.65.Ge Solutions of wave equations: bound states 
03.65.Nk Scattering theory 

21.60.Jz Nuclear Density Functional Theory and extensions 

1 Introduction 

As protons and neutrons arc indistinguishable in nuclei, Pauli principle must be taken into account in order to 
build nuclear wave functions. This has the well-known consequence of the appearance of exchange potentials 
in the one-body Schrodinger equation, in the optical potential of reaction theory or in the Hartree-Fock 
(HF) potential, issued from the variational principle applied to an antisymmetric wave function of independent 
particles [5]. Exchange potentials are non-local, i.e they are integral operators acting on the wave function 
at each point of space, so that Schrodinger equation written in coordinate space representation becomes an 
integro-differential equation. This type of equation cannot be handled by standard numerical methods used in 
ordinary differential equations, such as midpoint and Henrici schemes 13^. 

Many methods have been devised to deal with non-local potentials. While some of them consider the 
implementation of bound states only scattering states can be handled with the methods of Refs.[51 [SI [7], 
where the Schrodinger equation is transformed to an integral equation. Solving the Schrodinger equation 
represented in momentum space has also been considered (see Refs.[8j [9j), as it becomes an integral equation 
therein as well. However, in practice, momentum space has to be discretized [9], so that wave functions do 
not have proper asymptotic behavior when back-transformed to coordinate space. Moreover, consideration of 
Coulomb potential cannot be handled exactly, its Fourier-Bessel function being undefined due to its infinite 
range. In the FRESCO reaction code (see Ref. 10 for numerical methods employed therein), dealing with 
coupled-channel sets of equations involving non-local couplings, non-local parts in equations are replaced by 
source terms, converging iteratively to their exact values, which is the most straightforward method to transform 
a non-local equation into a local equation (we will mention it from now on as the source method). However, the 
source method usually converges very slowly, and sometimes even diverges. Indeed, even in the case of reaction 
theory, where local potential is dominant, Pade approximants have to be used in the FRESCO code in order to 
avoid instabilities ^ . 
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A very simple method has been introduced in Ref. [TT] , where the non-local potential is replaced by a state- 
dependent potential, the so-called trivially equivalent local potential (TELP), and where the problem has to 
be solved iteratively as well. This scheme is particularly interesting in HF framework, where iterative scheme 
is the only one available as HF equations are non-linear. As all considered equations become local, standard 
numerical methods for ordinary differential equations can be used, so that the implementation of both bound 
and scattering states poses no problem. Moreover, it is very stable in practice, contrary to the source method, 
so that it is widely used in nuclear theory [141 [151 [16] . However, as the TELP method involves a division by the 
considered wave functions, it contains poles, which are difficult to handle numerically. The solution proposed 
in Ref.|llj consists in an interpolation of the potential near its poles, replacing the potential by a straight line 
connecting two points before and after the pole, the distance between the two points going to zero iteration after 
iteration. However, this procedure is difficult to apply, as there is no precise criterion to determine how fast 
the mentioned distance has to decrease on the one hand, and, on the other hand, potentials become very large 
close to poles, which can generate numerical inaccuracies. It would then be interesting to have the advantages 
of both source and TELP methods, i.e. fast convergence of the iterative scheme and absence of divergences in 
potentials. 



2 Combination of source and TELP methods 



2.1 Local equivalent equation 

For simplicity, we will consider spherically symmetric potentials only, even though the method can be readily 
extended to more complicated situations, such as coupled-channel sets of equations. The Schrodinger equation 
then reads: 
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where is the charge of the target (proton) or is equal to zero (neutron), Cc is the Coulomb constant ~ 1.44 
MeV fm, m is the effective mass of the nucleon in MeV c~^ units, u{r) is the calculated wave function, of 
orbital momentum and energy e in MeV units, v{r) is the remaining local part of the potential in MeV units 
and w{r,r') its non-local part in MeV fm~-^ units. v{r) and w{r,r') are assumed to decrease quickly at large 
distance, so that u{r) becomes a linear combination of Coulomb wave functions for r > R, with R sufficiently 
large. 

We will now define the equivalent local equation used in the proposed method: 
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where Ubefir) is the normalized wave function obtained at previous iteration, vioc{r) is the state-dependent 
local potential used in our method, s{r) is a source term, F{r) is a smoothing function, detailed afterward, and 
Cbef is a normalization constant defined so that Ubefir) ~ Cbef f^'^^ for r ~ 0. Boundary conditions demanded 
for u{r) in Eq.([2|) are u{r) ^ r^~^^ for r ^ 0, and outgoing wave function condition at r — > -|-oo if one considers 
a bound or resonant (Gamow) state [12|, I13j. The boundary condition for u{r) at r ~ implies the presence 
of Cbef in Eq.dH), because, during integration of Eq.®, u{r) is not normalized whereas Ubef{r) is. u{r) will 
naturally be normalized at the end of the calculation. If u{r) is a scattering state, no boundary condition at 
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r +00 is required. One can easily check that Eq.Q is equivalent to Eq.(IT]) if Ubef{r) oc u{r), hence at 
convergence of the iterative process. 

2.2 Interest of the method 

The method embodied by Eqs. (|2|3|4|5p is very close to the TELP method of Ref. [llj . TELP method is indeed 
recovered if one arbitrarily sets F{r) = in Eqs. (|3|4p . The presence of F{r) is hence directly related to the 
zeroes of the function u{r), which are the cause of the divergences in the TELP method. F{r) consists in two 
factors in Eq.([3]). The first one is readily seen to be virtually zero except in the vicinity of the nodes of u{r), 
where it behaves like a Gaussian of maximal value equal to one, so that it cancels the divergences present in the 
TELP. The second factor originates from the fact that the TELP is well behaved at r '--^ 0, even though u{0) = 
[TT] . Thus, it is numerically more stable to have vioc{r) equal to the TELP close to r = 0. As the second factor 
is virtually equal to zero close to r = 0, whereas it is otherwise almost equal to one, it removes the action of the 
first factor at r = 0, while leaving it unchanged for the other nodes of u{r). The used decay constants equal to 
100 are rather arbitrary and were empirically determined in order to have stable calculations for a large set of 
nuclei. As a consequence, vioc{r) is always finite close to the nodes of u{r). Hence, the integration of Eq.(l2]) is 
always numerically stable. Moreover, as s(r) is non zero only in very small regions of considered radii, it does 
not hinder convergence as in the source method. Note that u'{r) does not enter Eq.(l2]), so that Numerov and 
Henrici methods can be applied to solve it. Consequently, both advantages of the source and TELP methods 
are present in the proposed numerical scheme. Another interesting feature is the possibility to obtain rapidly 
approximate energies of bound and narrow resonant states of Eq.Q. For that, one diagonalizes vioc{r) with a 
basis of harmonic oscillator (HO) states. Due to the overall smallness of s{r), eigenvalues of the HO matrix are 
very good starting energies, which are refined afterward with Newton method in order to determine the exact 
eigen-energies of Eq. ^ . 

3 Example: HF potential of ^^O 

3.1 Motivation 

In order to illustrate the effectiveness of the technique described in Sec.©, we will consider the evaluation of the 
HF potential of ^^O generated by the realistic interaction N'^LO [T7], renormalized within the low-momentum 
interaction framework [18] . The maximal momentum in two-body relative space used for the latter is A = 1.9 
fm~^. Note that we do not aim at describing properly ^^O properties at HF level, but simply at illustrating the 
proposed numerical scheme. HF equations are solved iteratively using linear mixing method for HF potential 
and a Woods-Saxon potential was used as starting point for the HF iterative process. As we consider a rather 
small nucleus, a very large potential mixing of 80 % could be used, which resulted in a very quick convergence 
to the HF solution in 32 iterations. If pairing interaction is used, in the context of Hartree-Fock-Bogolyubov 
(HFB), linear mixing might be insufficient, so that the more powerful modified Broyden method should be used 
(see Ref. |19 for recent application to HFB formalism and comparison to linear mixing method). Moreover, 
as the N^^LO interaction is decomposed in a HO basis [20], the non-local part of the HF potential is a sum of 
separable functions of the form /(r) g{r'). As convergence of shell model energies and eigenvectors is very quick 
with the number of HO basis states, which was shown in Rcf. 20 , 9 HO states per partial wave are used in the 
decomposition of the N'^LO interaction in the present work. The used HO parameter was 6 = 2 fm. Hence, 
the integrations involving w(r, r') in Eqs. pi4p for each r are replaced by sums over a few number of HO states, 
rendering wave functions determination very fast even though non-local operators are used. Note nevertheless 
that potential separability is not demanded in our method, this property being used only to improve efficiency 
of calculations. 

3.2 Results 

Energies of occupied single particle states and of unoccupied lsi/2 states are shown in Tab.lH]), where one can 
see that the proton lsi/2 state is resonant. We will now concentrate on the proton and neutron Si/2 states, 
as they contain nodes in the nuclear region, which would induce divergences in TELP. Besides bound/resonant 
Si/2 states, 30 scattering si/2 proton and neutron states have been calculated (see Tab.© for the value of their 
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Table 1: Single particle energies of occupied Osi/27 OP3/2 and Opi/2, and unoccupied lsi/2 proton (Ep) and 
neutron (£"„) HF states of ^^O. They are given in MeV units. For the resonant lsi/2 proton state, width is 
written between parentheses in keV units after its energy value, 
state Ep En 

Osi/2 -64.657 -69.444 

Op3/2 -29.221 -33.821 

Opi/2 -19.753 -24.065 

lsi/2 0.519 (22.928) -2.826 



Table 2: Linear momenta of scattering si/2 proton (fcp) and neutron (fc„) HF states in fm ^. Used principal 
quantum number is arbitrary. 

S1j9i1jG /bp 

2si/2 0.105-10.00469 0.105 
3si/2 0.123-10.0231 0.123 
4si/2 0.15-10.05 0.15 
5si/2 0.177-10.0769 0.177 
6si/2 0.195-10.0953 0.195 
7si/2 0.214-10.0953 0.214 
8si/2 0.269-10.0769 0.269 



9si/2 0.35-10.05 0.35 

lOsi/2 0.431-10.0231 0.431 

llsi/2 0.486-10.00469 0.486 

12si/2 0.523 0.523 

13si/2 0.616 0.616 

14si/2 0.75 0.75 

15si/2 0.885 0.885 

I6S1/2 0.977 0.977 



linear momenta). As the proton lsi/2 state is resonant, proton si/2 scattering states have been chosen to belong 
to a complex contour in /c-space, as would be the case in a Gamow Shell Model calculation [T51 [TBI \Tn \n[ I23j . 
where Berggren bases of complex energy-states consist in bound states, resonant states, and contours of complex 
scattering states enclosing resonances [2^. Potentials vioc{r), sources s(r) and wave functions u{r) are illustrated 
in Fig.([T]) for bound neutron lsi/2 state and resonant proton lsi/2 state, and in Fig. ([2]) for two si/2 proton and 
neutron scattering states of linear momentum k = 0.977 fm^^. One can see that potential viocij) and source 
s(r) vary very much close to the nodes of associated wave functions, but their maximal values in modulus remain 
sufHclently small not to generate numerical inaccuracies. Moreover, these variations partially cancel in Eq.Q, 
as v{r) and w{r, r') potentials are smooth in Eq.(IT]). In order to check the accuracy of obtained wave functions, 
we have calculated the overlaps between the considered Si/2 states, in both proton and neutron cases. As one 
deals with unbound states, complex scaling method ^22J is used to calculate radial integrals, which diverge on 
the real axis. Average and maximal norms (used norm is max(|5R(z)|, |3(z)|)) of overlaps are shown in Tab.Q. 
As average values of overlaps are of the order of 10^^, calculated states are numerically orthogonal, which shows 
that the method described in Sec. ([2]) is very precise. 

4 Conclusion 

Numerical methods allowing to solve Schrodinger equations with non-local potentials are usually complicated 
or restricted to a given class of subproblems, such as the consideration of bound states uniquely. Source and 
TELP methods, on the contrary, are very simple to code, but at the price of slow convergence and reduced 
stability for the former, and of the appearance of potential divergences in the latter. The proposed method. 
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Figure 1: Wave functions (up), potentials (middle) and sources (bottom) of proton and neutron lsi/2 HF states. 
Real part (Re) is depicted as full line and imaginary part (Im) as dotted line. Imaginary part has been multiplied 
by a factor of 1000 in the middle/left and bottom/left quadrants to be visible on the figure. Imaginary parts 
are present for the proton case as proton HF state is resonant, whereas the bound neutron state is real, so that 
no imaginary part occurs therein. See Sec. (|2.ip for definitions of potentials and sources. 
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Table 3: Maximal and average norms of the overlaps between the considered Si/2 states, i.e. bound/resonant 
Osi/2 and lsi/2 states, and the scattering states described in Tab. ([2]). Overlaps are calculated with the complex 
scaling method [22]. 





Maximal 


Average 


Proton 


1.526 xlO"* 


2.207 xlO-9 


Neutron 


1.287 xlO-* 


1.307 xlO-9 



by combining advantages of both source and TELP methods, allows to solve non-local Schrodinger equation 
very quickly and precisely, which has been shown with the example of HF calculation of ^^O with finite-range 
realistic interaction N'^LO. Due to its simplicity and efficiency, it is a very interesting method to deal with the 
integration of non-local Schrodinger equation. 
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